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We construct an algorithm to simulate imaginary time evolution of translationally invariant spin 
systems with local interactions on an infinite, symmetric tree. We describe the state by symmetric 
infinite- Tree Tensor Network (iTTN) and use translation-invariant operators for the updates at 
each time step. The contraction of this tree tensor network can be computed efficiently by recursion 
without approximations and one can then truncate all the iTTN tensors at the same time. The 
translational symmetry is preserved at each time step that makes the algorithm very well conditioned 
and stable. The computational cost scales like 0(D q+1 ) with the bond dimension D and coordination 
number q, much favourable than that of the iTEBD on trees [D. Nagaj et al. Phys. Rev. B 77, 
214431 (2008)]. Studying the transverse-field Ising model on the Bethe lattice, the numerics indicate 
a ferromagnetic-paramagnetic phase transition, with a finite correlation length even at the transition 
point. 
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I. INTRODUCTION 

Since S. White's Density Matrix Renormalization Group method [l|, Q was established, many 
developments have been made in order to simulate interacting spin systems as they offer an efficient 
way to study the properties of their ground states and are applicable to fermionic systems as well, 
see @ for a review. The matrix product states description gives us an accurate approximation of 
the ground states of gapped Hamiltonians on a spin chain, but it can also be applied to detect phase 
transitions. Various numerical methods have been proposed to carry out these simulations and have 
achieved decent success. Their higher dimensional generalizations also appeared the PEPS @, 0] and 
its extension to infinite lattices, the iPEPS @ and @, tree tensor networks [HJ, [HI, MERA pjf and 
other approaches are promising, but their extensive use is hampered by their large (but nevertheless 
polynomial) computational costs. 

In one branch of these methods, that uses MPS description in one spatial dimension and more 
generally tensor networks in higher dimensions, one describes the state of a spin system in a lattice 
by placing tensors to each vertex. These tensors has one 'physical' index corresponding to the local 
spin basis vectors, and 'virtual indices' for each bond coming out from the vertex. The dimension of 
this index, the 'bond dimension' is a parameter. The quantum mechanical amplitude of a given spin 
configuration is obtained by fixing the physical indices and contracting all the virtual ones. In one 
dimension, the computation of the expectation of a local operator O i.e. (ip\0\ip) can be interpreted 
as the contraction of the chain since we have to sum over all the virtual indices. This can be computed 
exactly for a fixed bond dimension in ID, but for higher dimensional lattices this does not hold any 
more. To see this, consider a subset of the lattice and contract all the indices belonging to the 'inner' 
bonds. Clearly the resulting tensor has a number of indices proportional to the length of the boundary 
of the subset. Several approximative techniques have been developed (see e.g. [fj, [HJ]) to overcome 
this and carry out the contraction in polynomial time. 

In this paper, we propose an algorithm to simulate imaginary time evolution of a translation- 
invariant Hamiltonian with local interactions on the Bethe lattice. We essentially generalize the 
matrix product operator (MPO, [l5j]) description for higher dimensional lattices and use the fact that 
a translationally invariant infinite tree tensor network on the Bethe lattice can be contracted efficiently 
by recursion. This follows directly from the geometry of the tree and is essentially the generalization 
of the power method for finding the leading eigenvector of a matrix, but instead of the matrix we 
have a tensor with q indices. This allows us to update all the iTTN tensors at once in every time 
step. Thus, all the symmetry is preserved during the whole evolution which makes the algorithm 
numerically very well conditioned and stable. 

Besides its symmetry preserving and numerical stability, the biggest improvement over the previ- 
ously established iTEBD method for trees [3) is the considerably lower computational cost. In [T^ |. 
this scaled like 0(D S ) for a Bethe lattice with coordination number q = 3, whereas for the present 
approach it is 0(D 4 ). 

Then we test our new method by simulating the ground state of the transverse field Ising model 
defined on the infinite tree. Plotting the transverse and longitudinal magnetizations we find a phase 
diagram qualitatively similar to that of the Ising model on a spin chain with distinct ferromagnetic 
and paramagnetic phases. The energy and its first derivative with respect to the external field are 
continuous, while the second second derivative has a jump. Surprisingly, the numerics indicate a finite 
correlation length for any external field, a striking difference from the one-dimensional counterpart. 
Remarkably, we find a very good agreement with previous findings by the iTEBD method [l4| and 
path integral Monte Carlo fig . 



II. CONTRACTION OF TREE TENSOR NETWORKS 
A. Tensors 

Consider a tensor A. The number of indices is sometimes called the order or degree of A. Although 
misleading, it is sometimes even called the rank, for example A a p 1 is a 3rd order tensor (or sometimes 
a 3-way tensor), a 2nd order tensor is a matrix, a first order tensor is a vector. 
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We will say that A is fully symmetric, if A aia2 „ an = A a , 1)a , 2 , . a . } for any permutation n of the 
set {l,2...n}. 

Multiplications of tensors make sense if we specify the indices to be contracted. For instance we can 
define the j-mode multiplication of the tensor A € ^{hx..Jj...xi n } by the vector v € C /j producing a 
tensor A Xj v e CU 1 *" 1 *- 1 ' 1 ^ 1 "* 1 *} as follows 

(A X j ujctj ,,,<Xj—\ ,(Xj+x ...<x n ■ — ^ ^ Act 1 ...a : j...a n 'Uctj 

where we also introduced the symbol Xj as in Q. Similarly one can multiply A by the matrix 

M <= C^ XJ i, as 

[A Xj Af) Ql ... Qj .... Qri := j4 ai ... a ' a n M ajCl / . 
a ', 

We will also make use of the outer product, i.e. one can make a 3rd order tensor A out of vectors 
u, v, w as 

A a /3-y = (u O V O w) a /3~ f := UaVpWj. 

B. Contraction, expectation values and correlations 

Consider the Bethe lattice (Fig. [IJ with coordination number q and place fully symmetric TTN 
tensors A s to each node. This describes the translation-invariant state ip; from each node, the tree 
looks the same. 




FIG. 1. Three layers of the Bethe lattice with coordination number q — 3 

Observe that the computation of the norm square (V'lV'} on trees essentially reduces to the compu- 
tation of the leading vector r of the transfer matrix E = J2 S A s ® A s , and it can be done efficiently 
by recursion. This statement is well known and easy to see in ID, but similarly simple also for trees. 
Without loss of generality take q = 3 and use tensor notations for simplicity. Begin with the 3-way, 
fully symmetric, d x d x d tensor E and a random vector A ^ of dimension d. For n — 0, 1, 2... define 
the series 



r («+i) =£ Xl r (n) x 2 r (n) . 



(1) 
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The normalized series converges to the so-called leading vector of E, i.e. r( n )/ll r ^ll r ' gi ven 
( r l r o) 7^ so the leading vector r is a fix point of the recursion (JT|) 

r = E Xi r X2 r. 

Notice that this is just the extension of the power iteration, which is originally an algorithm for finding 
the dominant eigenvector of a matrix. Its generalization for the case where E is not symmetric and 
has different leading vectors for each mode is straightforward. Notice also, that for a generic 3rd order 
tensor T, its leading vectors u,v,w give the best rank-one decomposition [ij], i.e. minimize 



IT - Ait 



O V o w\ 



and A = rx 1 MX 2 iix 3 iij, assuming ||u|| = ||u|| = ||w|| = 1. 
Using the leading vector r one gets 



(i>\i>) = (r\r) =y^r 



2 

r 



The expectation of a local operator O in the state ip reads 

(r\E \r) 



iO) = 



(r\E\r) 



with Eo = J2 S S '( S \0\ S ')A S ® A s ' , and the correlation between two sites of distance n reduces to the 
computation of the second largest eigenvalue of the matrix 

M = E x i r. 

This is in accordance with the findings of [14( but written in a much more compact form. 



C. Computational costs 

Clearly, the direct computational cost of each iteration step scales as 0(d q ) for arbitrary q € N. As 
for the convergence, we can state the followings. If the fully symmetric E can be decomposed such 
that 

E = XiVi o Vi o Vi 

i 

where (vi\vj) = for i ^ j and |Ai| > |Aa| > | A3 1 > ... then the convergence is surely geometric with 
ratio I A2/A1 1 . This corresponds to the matrix case when the matrix is diagonalizable. In any case, 
given a concrete E, the convergence should be checked by doing the numerics. Let us mention that 
since IA2/A1I is related to the correlation length of the system, one can expect an exponentially fast 
convergence only if the correlation length is finite. Otherwise, if the correlation length diverges as the 
function of the bond dimension, the algorithm slows down. 

However, we can exploit the symmetries of E to get a much better scaling. To this end, simply 
exploit the structure of E = J2 S A s ® A 3 . Rewrite ((T|) explicitly 

1 

(n+l) _ V V 4 s I s r {n) v (n) (?) 

where all the Greek indices run from 1 to D. Now break the sum into four parts. First let us 
sum over a, then over /3. When regarding r aa i as a matrix, we will write R, i.e. define the tensor 
B s = A s x 1 R( n > , explicitly written 

D 

°q'^7 ■— / , A a8-y r aa'- 
a=l 
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This requires 0(D 4 ) steps because we fix s, a', f3, 7 and sum over a. In the same way perform 

D 

C s — K s r {n) 

^a'P'-y ■— / , a u'8n' 88" 
/3=1 

Finally we arrive at 

1 D 
X! ^2 C a'8>~/ As a >8>y 

s=0 a'8' = i 

again at a cost of 0{D A ). Performing the recursion as described above requires 0(D A ) computational 
steps. Its easy to see that the idea works for arbitrary q and its cost scales like 0(D q+1 ). 

Note that, the above argument still holds for the computation of the expectation value (O) . 

In order to compute the correlation length one has to find the second largest eigenvalue of the 
D 2 x D 2 matrix M at a cost of 0(D 4 ) steps. 

D. Tensor updates 

Let us consider a time evolution induced by a translationally invariant Hamiltonian on the tree lat- 
tice. If one starts with a translation-invariant state, then the symmetry of the Hamiltonian guarantees 
that it is never broken in time. Assume we have a TTN description with bond dimension D for this 
translationa- invariant state tp. Obviously each site has the same TTN tensors A s . Now we update 
these tensors in every time step and obtain the updated state ip' with new TTN tensors and bond 
dimension 2D. Then we need to truncate them in an optimal way, i.e. we should find TTN tensors 
A s with bond dimension D describing a state if) which minimizes the distance from tp' , i.e. for which 

livz-^ii 

is minimal. 

In order to carry out the optimal cut express the norm square by the leading vector r' of the updated 
E' as 

2D 

WW) = E' X! r' x 2 r' x 3 r' = r'xr'= ^ r'^r'^, = {7\r'). 

77'=1 

One can look at r^, as a symmetric 2D x 2D, complex matrix which we will denote by R' = R'^ . 
Perform its singular value decomposition: R' = UAW and introduce the D x 2D matrix U consisting 
of the first D rows of U and A is a D x D diagonal matrix made up of the D largest singular values. 
Define the self adjoint projection P — U^U, indeed P = P^ = P 2 since UW = 1 D . We can project 
to the subspace of the highest D singular values as A s x 1 P x 2 P x 3 P, implying that the appropriate 
cut is then given by the transformation 

A s — >A s = A s x 1 Ux 2 Ux 3 U. 

These new tensors will have the right size D x D x D and describe the state ip. This is nothing 
else than the extension of the ID algorithm to the tree geometry. Recall that in one dimension the 
truncation step looks like [l5[ 

A s — > I s = UA S W = A s xi U x 2 U. 

As always, the lower the discarded weight relatively, the higher the accuracy of the approximation. If 
the former is not small enough, one should increase the bond dimension. 

Translation-invariant iTEBD on the Bethe lattice 



Jn+l) 
77' 
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FIG. 2. The norm condition Q in diagram 

We also show that the one dimensional iTEBD algorithm fig which uses the so-called canonical 
representation of infinite-MPS (iMPS, see e.g. the review Q) can also be extended to the Bethe 
lattice. Note that here one updates all the tensors at once therefore translation invariance is kept 
at each step, given that the infinite-MPO (iMPO) and the starting state have been so. Hence this 
update method differs from the one in [l4| where instead of iMPO, the interaction between only two 
adjacent sites has been considered at each step, leading to the breaking of the translational symmetry 
during the process. 

Canonical iTTN 

This representation contains the tensors T, placing to the sites and the explicit Schmidt coefficients 
A (sitting on the bonds) that one gets if cutting the tree between two neighbouring sites [3| 

i^>= n e a ^ n e r»& lt j..>i-i>i..) 

fcGbonds a k iEsites Si ,oti . . . 

each index a/ appears in two T tensors and one A. The normalization conditions for the canonical 
iTTN on the Bethe lattice look like 

E A « = 1 > ( 3 ) 



2-^2-^2-^ 1 a k ,ai,a' m A a k eti 1 a k ,a u a m ~ «m,V V 1 ) 
« a k ai 

this latter relation becomes much more unambiguous if represented in diagram, see Figj^l There 
are two conditions like this if q = 3, interchanging the legs of the tensor T, i.e. we have two layers 
due to the tensor product, and we leave two open legs (indices) in one direction and close the other 
legs by placing A 2 and sum over the corresponding bonds. The directional symmetry manifests itself 
in a fully symmetric V s . Now applying the iMPO to this iTTN results in higher bond dimension 
and the violation of these norm conditions: The new leading vector will be different from the identity. 
According to the procedure in (l6j , we have to bring this back to the canonical form (orthogonalization) 
and then truncate all bond indices by retaining the D highest Schmidt coefficients. We can naturally 
extend this by performing the same operations on each leg as in the one dimensional case, noting that 
we can define the dominant eigenvector by (JTJ) and in case of directional symmetry, this will be the 
same for all legs. After computing this leading vector R (which is interpreted as a matrix because of 
the two layers), we perform the same manipulations as in Fig. 3 in [la ], but now we have open legs 
in one direction and 'closed legs' in the other q — 1 directions. 
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III. IMAGINARY TIME EVOLUTION 

Consider the Hamiltonian on the Bethc tree G = (V, E) with coordination number q 

As shown in THJ one can derive an elegant MPO representation for e~ eH with H being the Ising 
Hamiltonian on a one-dimensional infinite chain. 

exp(-eff)= Tr{C kl (£)C k2 (e)...)Z kl <g> Z k \.. (6) 



with 2 dimensional MPO matrices C k (e) given by 

<%> = I>«f = ( c t\iL) 



r<\t n D d t / Vcoshesinhe \ 

C ( £ ) = ^S iffil B i - ^ Vcosh£sinhe ■ 

i 

And the spin operators are A = e Sa , A 1 = a x , S = hAt/2 for the second order Suzuki- Trotter 
formula [l5(. We would like to determine the local 'transfer tensors' Q s Jp 1 that makes the TTN tensors 
evolve as 

s' 

or using tensor notation 

A s — > J2A S ' ®Q s ' s . 

s' 

We can express the MPO matrices by the sum of outer product of vectors. First define the vectors 



v = (^cosh(JAt), y/sinh(JAt)) 

and 

u = (cxp(/iAi/2), exp(-/iAt/2)). 
Then observe that in ID, the MPO-s has the compact form 

= v a v p l{k = a + 13} 

a,/3 corresponding to the x direction, and 1{^4} being the indicator function of the event A. 
The translationally invariant local transfer tensors read 

Qa'p = V a Vf}U s U s <l{s + s' = a + /?}. 
analogously, if there are three spatial indices, the MPO-s has the form 

C k M = v a v,3V 7 , l{a + P + 7 = k} 
and the translationally invariant local transfer tensors read 

Qa/3j = VaVf3VyU s U s/ l{s + s' = a + /? + 7}. 

this can be written as a sum of two rank-one tensors as can be verified 

Q = uouovovov J t-iiouovovov 

using tensor notation and u = a z u, v — a z v. 

Note that the above derivation for the local interactions like a x ® a x with a transverse-field can be 
repeated for a y ® a y or a z ® a z . 
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IV. RESULTS AND DISCUSSION 




FIG. 3. Longitudinal and transverse magnetization of the Ising model on the q = 3 Bethe lattice. 
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FIG. 4. Longitudinal magnetization of the Ising model on the q — 3 and q — 4 (subplot) Bethe lattice around 
the transition point. Notice that even for small bond dimensions, the results are barely different. 

We set J = 1 making interactions ferromagnetic and vary h. 

Using the present algorithm I found very good agreement with the results reported earlier. Namely, 
the transverse and parallel magnetization as the function of the magnetic field behave similarly to 
that of the Ising model on an infinite line, but the critical point is shifted, and estimated to be in the 
interval 2.23 < h/J < 2.25. 
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FIG. 5. Correlation length around the transition point for different bond dimensions, q — 3 



One can also check whether the longitudinal magnetization obeys the scaling law 

(a x ) oc (h c - h) 13 

in the ferromagnetic phase, near to the transition point. Fitting a line in the log-log plot, we obtain 
8 w 0.46. In [14j they reported 8 = 0.41, while the so-called 'cavity method' prediction 8 = 1/2 
equals the mean-field result, see jl8j ]. 








FIG. 6. Log-log plot of the parallel magnetization vs h — h c 

Surprisingly, the correlation length £ at the transition point does not seem to diverge as D — ¥ oo, 
since the second largest eigenvalue A2 of B never exceeds 1/2, in complete agreement with the findings 
of [l3|. This is a strong numerical evidence of a finite correlation length even in a presence of a 
(supposed) phase transition. Also, plotting the first derivative of the energy per site with respect 
to the magnetic field, we see a non analytic breaking-point at the suspected critical field, while the 
plot suggest a discontinuous second derivative with a finite jump at the same point, see Fig. [7] 
Qualitatively the ground state energy behaves very similarly to that of the exactly solvable quantum 
Curie- Weiss model [l8j . 

Note that in [lij they used a different parametrization of the Hamiltonian with the parameter s for 
which h/J — 3(1 — s)/s. 

I also considered the Bethe lattice with coordination number q = 4 and found similar phase dia- 
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grams, but with transition point at h/ J ~ 3.3 and the results suggest that A2 < 1/3. 
Previous results and theoretical considerations 

In [l7| the authors extended the cavity method to quantum spin-1/2 models making it applicable to 



Bethe lattices, too. In 18], by using Monte Carlo simulations, it is applied to the present model. They 



computed the magnetization for different temperatures, and their zero temperature extrapolations 
(e.g. for the ground state) agree very well with the present findings and with [3] ■ 

We would like to draw the reader's attention to the last section of [l4j], where the two point cor- 
relation is calculated. Given a translation-invariant Hamiltonian with nearest neighbour interactions 
on the Bethe lattice, it is shown that the two point correlations always fall off exponentially, and the 
correlation length is upper bounded by l/log(q— 1), confirming our numerical findings. According to 
this brief argument, it is a model-independent and computational method independent consequence of 
assuming that the translational independent ground state is the stable limit of a sequence of ground 
states as the size of the tree grows. 

Notice that the argument also holds for the classical Ising model defined on the Bethe lattice which 
is exactly solvable and exhibits phase transition (l9j . 



Note also, that in [18j, the correlation length was not computed, but (mistakenly) supposed to be 



infinite at the critical point. 
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FIG. 7. Ground state energy and its first and second derivatives with respect to h for q - 
results for D = 8 and D = 16 are indistinguishable in the figure. 



3. Notice that the 



Errors 

Concerning the errors one can distinguish two sources: using the Trotter-Suzuki formula and the 

truncation. The former has the order of N(At) 3 = t 3 /N 2 , as the second order Trotter formula is used, 

ft t t \ N 

i.e. e~ tH ps I e _ 2iv^e _ "« ffx e~2N Hl J . So we need t — > oo to obtain the ground state properties and 

TV — > oo together with At = t/N — > and N(At) 3 = t 3 /N 2 — > to make the Trotter error go to zero. 
This can be achieved e.g. by setting N(t) = ct 2 for some constant c and then increasing t until the 
results do not change significantly. 

In Fig. [S]we plot the longitudinal magnetization m x — (<r x ) computed as a function of time. Here 
we fix the maximum time T and run the algorithm five times with different number of time steps N . 
One can see that higher the N the slower the convergence to the limiting value m x (T,N), but as a 
function of N (or At = T/N) the convergence to the limit m x (T) = limjv^oo m x (T, N) is very fast. 
Then increasing the time m x (T) — > m x , in the present case we can choose N(T) — 2T 2 . In general, 
if the system is gapless, this convergence is exponential in T . 
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Note that one has to be cautious about the convergence of different quantities, usually the conver- 
gence of the energy is much faster than that of the magnetization, so it is not sufficient to check only 
the convergence of the energy. 




FIG. 8. Typical behaviour of a computed physical quantity during the time evolution. The maximum time 
T is fixed, and the time steps N are increased linearly. Starting from the same initial state, we get different 
curves due to the Trotter errors. 



V. SUMMARY 



In this paper we have derived an algorithm to simulate the imaginary time evolution induced by 
certain translationally invariant local Hamiltonians on the Bethe lattice. During the time evolution 
we have used a translationally invariant iTTN description of the evolving state and have given the 
optimal updates which retain the symmetry at every step. The presented method has much lower 
computational costs than that of the slightly modificated version of the iTEBD method for trees, 
reported in [T3]. That scaled like 0(D 8 ), and the translational symmetry was broken during the 
algorithm. Investigated specifically the Ising model on trees, we have found very good agreement with 
previous numerical results, namely a phase transition between a ferromagnetic and a paramagnetic 
phase but with finite correlation length even at the transition point. At this critical point, the second 
derivative of the energy with respect to the magnetic field has a finite jump, while the energy and 
its first derivative is continuous. Considering the results for q = 2,3,4 suggests that increasing the 
coordination number q, at fixed J, the critical magnetic field increases, and the maximal correlation 
length decreases. 
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